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■ Abstract 

a ■ 
o 

We investigate the Dirac time-dependent variational method using a Gaussian 
^ ' trial functional for an infinite one dimensional system of Bosons interacting 

through a repulsive contact interaction. The method produces a set of non- 
linear time dependent equations for the variational parameters. By solving 
the static equations we have calculated the ground state energy per particle. 
We have also considered small oscillations about equilibrium and obtain mode 
equations which lead us to a gapless dispersion relation. The existence of an 
exact numerical solution for the ground state energy and excitations obtained 
by Lieb allow us to compare with the Gaussian results. We can also, as the 
system becomes less dilute, see the improvement of the results as compared 
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with the Bogohubov scheme. 
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I. INTRODUCTION 



Recently Bose- Einstein condensation in atomic traps was achieved |T|- [0] in a dilute 
regime {a^n « 1). This is in contrast with the Helium 4 regime where a much higher 
density leads to {a?n ^1). For the high density regime there is a large "depletion" 90%) 
^ where as for the dilute regime it is very small (~ 1%) [^. 

Theoretically these new experiments have been described very successfully with mean 
field theory without quantum fluctuations using the Gross Pitaevskii equation |^ or the 
Bogoliubov scheme The Gaussian self-consistent approximation presented here should 
be useful at in intermediate regime. Unfortunately due to three body recombinations there 
seems to be a limit for increasing the number of particles in the system. Recently it 
has been pointed out |^ that using a strong magnetic field it would be possible to make 
the system effectively less dilute by drastically changing the scattering length. This makes 
the comparison between self-consistent results and dilute theories very important. In this 
context the one dimensional delta function case can produce some insight because a contact 
interaction can be used in the self-consistent theory in contrast with the 3 dimensional case 



IP - ||T0|- The existence of an exact solution for the ground state energy and for particle 
and hole excitations makes the comparison very intersting, provide that we understand how 
the separated particle and hole excitations of the exact solution [^| are connected to the 
particle-hole excitations given by the approximate methods. 

The objective of this paper is to exhibit the most general way of obtaining time dependent 



equations of motion in the Gaussian approximation [|r^. This will lead to the so called 
generalized RPA, when one examines infinitesimal oscillations about equilibrium. The static 
solution in the uniform case can be obtained using several other methods ||12[ - leading to 
a gap in the quasi-boson energy. We show here that the time dependent RPA equations lead 
to a gapless mode. In fact this must happen because particle number conservation symmetry 
is broken in the static solution, so the zero gap is exactly the associated Goldstone mode. 



This discussion can be seen as an alternative to the functional derivative [|T3|, method 
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in the Girardeau- Arnowitt approximation. 

The Bogohubov scheme, for a dilute or weak interacting system, can be obtained by 
a particular truncation of the Gaussian results. So that we can compare the Gaussian 
variational results , the dilute Bogoliubov scheme and the exact solution, for the particular 
case under discussion here. 

The structure of this paper is as follows. Section 2 reviews the time-dependent variational 
principle and the canonical nature of the equations of motion arising from it. In section 3 
we specialize to the one dimensional uniform case end examine the ground state energy and 
the excitations for both the approximate methods and the exact solution. Section 4 contains 
our numerical solutions and conclusions 



II. GENERAL FORMALISM 



In this section we shall review some of the results of the Time- dependent Variational 



Principle |TT[], [|T6[ and show how it can be implemented in the non-relativistic case. First 
we define an effective action functional for the time-dependent quantum system 



S = jL{t)dt = J dt{^,t\{tdt- H)\^,t), 



(1) 



where l^', t) is the quantum state of the system and H is the Hamiltonian of the theory. For 
a system of non-relativistic interacting Bosons we have [we use the notation : /x = / d^x] 



H 



V^(x)t/i(x, y)^(y) + \l ^(y)t^(x)tV(x - y)^(x)V^(y) 



(2) 



where the one body Hamiltonian /i(x, y) may include a one body external potential. The 
creation and destruction operators ip"^ and ifj can be written in the form 



^(x) 



x/2 



(x) + 27r(x) 



(3) 



V'(x)t 
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X — 27r X 
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where 0(x) is the field operator and 7r(x) the canonical field momemtum 

We can obtain the time dependent Schrodinger equation by requiring that S is stationary, 
supplemented by appropriate boundary conditions, under the most general variation of |^, i). 
The variational scheme is implemented by chosing a trial wave functional describing the 
system. Working in the functional Shrodinger picture we replace the abstract state 1^1/, t) by 
a wave functional of the field 0'(x) 

(4) 

The action of the operators 0(x) and the canonical momentum 7r(x) are realized respectively 

by 



0(x)|*,i)^(/)'(x)*[0',i] 



(5) 



The mean value of any operator is calculated by the functional integral 

{^,t\0\^if,t) = (6) 

where "if is normalized to unity The Gaussian approximation consists in taking a Gaussian 
trial wave functional in its most general parametrization 



^^f[^',t] = Nexpl- f 50'(x,t) 



-iE(x,y,t) 



+ i£7r(x,i)5(/)'(x,i)|, (7) 

with 50'(x, = 0'(x) — 0(x, f). Due to the fact that the Hamiltonian commutes with the 
number of particles N = J^ip'^{x.)ip{x.) i.e. 

[H,N]^0 (8) 

we can actually define a more general trial functional 
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|^',t) = e-^^^W|^,t), (9) 

where 6{t) is another variational parameter introduced because of this continuous symmetry. 
Thus our variational parameters are 0(x, t), 7r(x, t), 6(t),G{:K, y, t) and S(x, y, t) ,with G and 
S being real symmetric matrices. These quantities are related to the following mean- values: 

(vi/',t|0(x)|^',t) = 0(x,t) 

(vl>',t|7r(x)|vl/',t) =7r(x,t) (10) 
(vj/', t|0(x)0(y) 1^', t) = G(x, y, t) + 0(x, t)0(y, t) 

(^', t^T-l^', t) = / S(x, y, t)G(y, x, t) + / 7r(x, t)0(x, t) 

+J\f9(t) + total time derivatives. (11) 

We may ignore the total time derivatives because they do not contribute to the equations 
of motion. If now we write the action we will get 

S = Jdt (^^ 7r(x, t)0(x, t) + ^ ^ S(x, y, t)G{y, x, t) + Me{t) - 7i) , (12) 

where 

7^ = (13) 

and 

Af = {^',t\N\^',t). (14) 

From ( P^ ) we see that (A/", 6), (tt, 0) and (S, G) are canonical pairs. Because of the symmetry 
7i has no dependence on 6 and it follows that = and d{t) = constant = fi . We can now 
write the remaining Hamilton equations, 

6{n- i^Af) 



(x,t) 



57r(x, t) 



S(x,y,t) 



5S(x,y,t) ' 
(5(7^ - /iAT) 
"5G(x,y,t) ' 
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For convenience we introduce 

V'(x,t)^(V^(x)) 
so that the equations for and vr become 

2V^(x,t) = 

To obtain Ti — fiAf we have to compute 



(x, t) + i7r(x, t) 

71 ^ 



6{n - ^M) 

(5V'*(x,t) 



using (^) and we have 

7^-/iAr = jww[(t^\t] L'(x) 



50'(x)) (^^'^^^ 50'(y), 



^(x-y) 
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(16) 



(17) 



(18) 



(19) 



All the functional integrals can be easily computed using an aditional source term (Ap 
pendix A) leading to 

n-^^^^= { [/.(x, y) - ^^5{^ - y)] p(x, y, t) + \v{^ - y) |^(x, t) |Xy, t) 

+ i 4y y (x - y) [i?(y, X, t)i?(x, y, t) + i?(x, x, t)i?(y, y, t) + D*(x, y, t)D(x, y, t)] 



+ /x,ynx-y) 



iv^*(x, t)^(y, t)i?(x, y, t) + ^^*(y, t)V^(x, t)i?(y, x, t) + |^(x, t)\''R{y, y, t) 



-- 4y y (x - y) [^(x, t)^(y, t)D*(x, y, t) + ^*(x, t)^*(y, t)D(x, y, t)] 



(20) 



and 



p(x, y, t) = (V'^(x)V(y)) = V*(x, t)V(y, t) + ^(x, y, t) 



A(x,y,t) = -(V'(x)z^(y)) = -z^(x, t)z^(y, t) + D(x, y, t) 



(21) 



with 



i?(x,y,t) 



G-i(x,y,t) 



+ G(x,y,t)-5(x-y) 



+ 



2 / S(x,w,t)G'(w,z,t)S(z,y,t) 

J w,z 



^(x,y,t) = 2 



+1 f [G(x,z,t)S(z,y,t)-S(x,z,t)G(z,y,t)] 

J z 



(22) 



G(x,y,t) 



+ 2 / E(x,w,t)G'(w,z,t)S(z,y,t) 



; 1^ [S(x, z, t)G(z, y, t) + G(x, z, t)S(z, y, t)] 



because of ([T0|) . It is easy to check that in terms of R and D the mean value Ti — fij\f 
corresponds to the standard mean- field factorization |jl2|. We note that the density gets 
contributions from the condensate field ip as well as from the fluctuations (GjS). The con- 
tribution from is the condensate density. So that the term with 4 ip^s can be interpreted 
as the condensate self-interaction. The interaction of particles not in the condensate with 
the condensate is taken into account by the terms with two ip^s. Finally the self- interaction 
of the particles not in the condensate comes from the terms with no ip {RR and DD). 
We introduce the generalized potentials 



Wd(x, y, t) = 5(x - y) / p(z, z, t)V{y. - z) 

J z 

We(x, y, t) = p(x, y, t)V{-K-y)= Ul + iUi 
Wp(x, y, t) = A(x, y, t) V(x - y) ^ + ^ 



(23) 



where the notation emphasises real and imaginary parts of Up. We also define the matrices 



^(x, y, t) = /i(x, y) - /i + W^(x, y, t) + W,^(x, y, t) + ^^(x, y, t) 
5(x, y, t) = /i(x, y) - /i - W^(x, y, t) + W,^(x, y, t) + ^^(x, y, t) 
C(x, y, t) = h{x, y) - fi + We(x, y, t) + Ud{x, y, t). 



(24) 



From Eqs. (05) and (|T^ we obtain an abstract matrix form of the equations of motion 
t 

G = {A, {G, S}} - [Ul G} - [Ul G] (25) 
= G^-Upifj*, 



^G-'AG-' - 2SAE - I + {Ul S} - [Ul S] 
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where we have used the fact that (S,G) are symmetric matrices. These equations (pSj) 
are the nonhnear field equations for an arbitrary interaction V between the particles and 
contain any external potential through h. As an example the matrix product G^^AG^^ can 
be written in coordinate representation as 



/ G'\^,z,t)Aiz,w)G-\w,y,t). 

J z.w 



(26) 



The static equations can be obtained by setting the canonical momenta to zero, that is 
S(x,t,t) = 7r(x,t) = 0, G'(x,y,t) = 0(x,t) = 0. From (|TBp and <^ we than have 

"G-Hx,y) 



/?(x,y,0) = i?(x,y) 
D(x,y,0) = D(x,y) 
^(x,0) = ^(x) = 



1 
2 
1 

" 2 
fx) 



G(x,y)-5(x-y) 



-G(x,y) 



(27) 



V2 



So that for the static case we have to self consistently solve 

\ [ G-i(x,z)A(z,w)G-i(w,y)-fi(x,y)=0 

4 Jz.w 



(28) 



5(x, z)^/'(z) - 2^/'(x)^2(z)\/(x - z] 



0. 



using (|2T| ) -(p^ ) . We note that if we constraint G = 1/2 for the static solution this leads to 
R = D = and A = B so that equation (44) is to the usual non linear equation for the 
single quantity |^ obtained from a many body product wave function (permanent) for 
bosons. However the time dependent eqs. (|25| ) are more general, because our trial Gaussian 
is actually a coherent state with an indefinite number of particles. 



III. CONTACT INTERACTION IN THE ONE DIMENSIONAL UNIFORM CASE 



A. Gaussian Approximation 



We will specialize the results of the previous section for the one dimensional case and for 
a contact interaction so that 
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V{x -y) = \5{x - y). 



(29) 



Because of the existence of exact numerical solutions we will treat the uniform case so that 
the momentum representation is the natural choice where the quantities A, B and G can 
simultaneously be diagonalized and (j){k) = (f)6{k). The static equations ( PBD become [we use 
the notation : = ^ JZo dk] 



-G-\k)A{k) - B{k) = 



B{0) - 



0. 



(30) 
(31) 



The static version of A and B from (p4D can be written as 



.4(k) = e(k) - ^ + Up + 2U + X- 



B{k) = e(k) - /i - Wp + 2W + 



(32) 



where e{k) 



2m 



and the generalized potentials become 



U = Ud=U^ = X I R{k') 

Jk' 

Up = X [ D{k') 

J k' 



(33) 



So that we can write the solution for f|30f) as 



and for (|3TD we have 



G{k) = 



m 

B{k) 



(34) 



or 



(35) 



B{0) = X4>^ 



(36) 



using (|3^ with (^) we can express D and R as functions of A, B 
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D{k) 



R{k) 



G-\k) 



G-\k) 
4 



-G{k) 



+ G{k)-l 



lB{k) - A{k) 
4 ^A{k)B{k) 

1 j B{k) + A{k) 
^\2JAik)B{k) 



(37) 
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From (132| ) , (|33| ) and (Q) we see that A{0) and -B(O) must be positive so that if A < equation 
(P^) demands that -B(O) < which is inconsistent with the previous statement. So the only 
possible solution in this case is = 0. For A > the symmetry breaking solution 7^ , 
using (|SBD, gives us 



/i = A / mk') - D{k')] + h\ 

J k' ^ 



(38) 



Having solved for yU we can rewrite A and B in (|3^) and (p3D as 



A(A;) = e(A;) + 2A / = e(A;) + 2Aa 



5(A;) = e(A;) + A02 = e(A;) + 2A6. 



(39) 



On the other hand using ( pTf ) and ( pQ]) we can write a pair of non-linear equations for a and 
h 



X r [b — a] 

2 A' ^[e{k') + 2\a] [e{k') + 2\b] 



(40) 



1 



e{k) + X[b + a] 



2 [ ^[e{k') + 2Xa] [e{k') + 2Xb] 
where we have used the total density constraint 



2 Jk' 



(41) 



which actually becomes our scale. This non-linear set of equations can be solved for a given 
p, given a and b. Once we have calculated a and b we can compute the chemical potential 
through 
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/i = A [2p - a - 6] . (42) 

In the same fashion the ground state energy density {E/L = (H)/L) can be computed 
obtaining 

E/L = ^a^ - - Xab + Xp" + K, (43) 

where K, the contribution from the kinetic energy can also be computed in terms of a and 
h as 

K=^-l eik') I ^ eik') + X[a + b] _ 1 

l^[e(A:')+2Aa] [e(A;') + 2A6] J 



2 Jfc 



As an aside we remark that for a dilute system we can approximate the self-consistent 
equations for a and b by truncating them at a second iteration. A first iteration on ( ^Of ) 
takes a ^ 0, b ^ p which implies zero depletion and Up = 0. This leads us to a non-pairing 
theory (Gross Pitaevskii equation 0). Then the next iteration leads to 

Xp f 1 



2 Jk' ^e{k>Y + 2Ape(A;') 

(45) 

b^p-'-[[^m±^-l 
2 Jk' \ ^e{k'Y + 2Ape(A;') 

Then we can calculate R and D truncating the self-consistency and giving the same results 
as the Bogoliubov scheme. Physically this means neglecting the effect of the terms that take 
into account the self-intercation of the particles not in the condensate {DD and i?i?). This 
approximation is usually valid for dilute systems where these terms are not important. With 
this truncation the ground state energy can be easily computed giving 



E X 

where the dimensionless parameter 7 is 



(46) 
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Returnig to our discussion we determine the excitations through the RPA equations 
which can be found by expanding all quantities around their equilibrium value [|r^. Thus 
we write 



G{k, k', t) = G{k)5{k - k') + 5G{k, k', t) 

E{k,k\t) ^ 6E{k,k',t) (48) 
0(fc,t) = (f)5{k) + 5(l){k,t) 
Tr{k,t) 5Ti{k,i). 

Thus we have written G and S in the basis where the equilibrium G is diagonal and kept 
terms up to first order in small quantities, of course the diagonal basis is plane waves. 
It will be useful to introduce new momentum coordinates so that 

P = k-k' (49) 

1 = ^ (50) 

and 

6G{k,k') ^ 6G{P,q). (51) 

We will see that P and q can be interpreted as total and relative momenta respectively of a 
pair of quasi-bosons. We can then write the RPA equations in a form where P is diagonal 
and can be considered as a dummy variable 

SG{q,P,t) ^K{q,P)6J:{q,P,t) + CK{q,P)67c{P,t)+ f S^iq, q' , P)6J:iq' , P,t) 

Jq' 

- 6t{q, P, t) ^M(g, P)SG{q, P, t) + CM(g, P)50(P, t) + [ Suiq, g', P)SG{q', P, t) (52) 

Jq' 

50(P,t) ^7r{P,t)A{P)+ [ CK(g',P)(5E(g',P,t) 

Jq' 

-6-k{P,t) ^<f){P,t)B{P) + f CM{q',P)SG{q',P,t). 

Jq' 

We note that for a given value of P the (vr, 0) degree of freedom is coupled to the much 
more numerous degrees of freedom (E, G) which are labeled by q. Different q values among 
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(E, G) are also coupled. Introducing the notation f{q' + P/2) — /| and f{q — P/2) — /_, 
we find non-diagonal matrices in {q, q') 



>SK(q, q', P) = A [G+ + G_] \G'^ + G'_] +\[G+- GJ\ \G'^ - G'_ 



^M(q,q',P) = ^ + ^ 



and diagonal elements 



G-^G 



i-<-<-i 



1-1 n 



G^ G\ 



+ 



(53) 



SK(q,P) = 2[A+G'_ + >l_G'+] 



G;^GlM+ + Gi^G;M_ 
SM(q,.f j 5 



(54) 



Finally we see the coupling elements between (tt, 0) and (E, G) 



CK(q,P) = A0 [G+ + G_] 



CM(q,P) = ^0 
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(55) 



which vanish when the symmetry in (j) is conserved (0 = 0). As pointed out above the 
equations are diagonal in P so we can interpret it as the total momentum of a pair of quasi- 
bosons. Because (5G, 5E and 50, ^tt are canonical variables we may invert the definitions of 
momentum and coordinate. For convenience, we define column vectors 



e(g,P,i) 



'^5S(g,P,t)^ 



V 



57r(P, t) 



n{q,p,t)^ 



^ 6Giq,P,t)^ 



I 



\ 



(56) 



Then we can write a coupled oscillator Hamiltonian that corresponds to the RPA equations 
of motion in a suggestive matrix element form 



(57) 
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where the matrixes M ^ and K are the generahzations of oscillator mass and spring constant 



K 



5k + sk ck ' 



V 



+ Sm Cm 



We may separate the diagonal part of i^RPA so that 



V 



Cm B 



(58) 



(59) 



where 



1 /sk 0\ /5S 

2 \ Q A \5 



'5G* 



TT 



1 / S'k Cm \ / 5S \ 1 

2 I CM / V 57r / 2 



Introducing the trivial canonical transformation 



sm 5S 



5S 



Stt ^ V B 5tt 



we obtain a simpler form for the diagonal part 

1 / smSk \ / 

2 \ AB J \6n 
If we define Qi and ^2 



50 



B 



1 



sm \ / 
B^ 

Sm Cm \ / 
Cm 



1 0~ 



-{6G* 6<p*) 
2 \ 1 



(60) 



(61) 



(62) 



fi2(g,P) = V^K(g,P)sM(g,P). 



(63) 
(64) 



and use the definitions of sk and sm from (U) we get, after some algebra, the remarkable 
result 



fi2(g, P) = JA+B+ + V^-fi- = ^^i(A;) + Qi{k') 



(65) 
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so that fti{P) and 1^2(0', P) can be interpreted as the one and two free quasi-boson energies. 
We note that ^2(0,-?) = 2Qi{P/2), which means that at zero relative momentum ^2(^,0) 
corresponds to two quasi-bosons with momentum P/2. Thus the oscillations of the 6(f), 
67c pair can be interpreted as a quasi-boson mode while the oscillations of 6G, 5S can 
be interpreted as an interacting pair of these same quasi-bosons. When = 0, we get 
ck = Cm = and the one and two quasi-bosons systems are calculated independently. 
When 7^ we must rediagonalize so that our final modes will be mixtures of one and 
two quasi-bosons. The variable q represents the internal motion of the quasi-boson pair 
with interaction given by the quantities S. In general this is a scattering problem and we 
must search for the scattering amplitude at a given energy and P, where the asymptotic 
conditions are determined by (0). In addition the coupling of one and two quasi-bosons 
will always lead to a bound state which is a particular mixture of the one quasi-boson mode 
with a bound component of the two quasi-bosons. 

As we did for the static results it is straightforward to see that the truncation that gives 
the Bogoliubov results implies neglecting the coupling {Hint = 0) and will lead to the usual 
result 

n\{P) = ^e2(P) + 2Ape(P). (66) 

We can see that the Bogoliubov excitations consider no interaction between the quasi-bosons. 
Note that because the Bogoliubov does not take into account the self-interaction of the 
particles not in the condensate we have k' = which means P = k. 

Returning to the dispersion relation for the bound mode one can finally eliminate the 
two quasi-boson components. To see this we try oscillatory solutions for (|52D such as 

e(t) = Ge*™ 

(67) 

n(t) = Ue'^\ 

and equations (^) can be written in a compact form 
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(68) 



that is 



V 



-Q 
B 



—Xx(f) 


^ 6G 







-A 


S(j) 







^ \{xX + rR) ^ 
-XyY - 2XzZ 

X(t)X 
-X(j){Y + 2Z) 



(69) 



where for simphcity we have used 



r{q, P)^G+- (?_ 



x{q,P) = G+ + G. 



y{q,P) 



1 + 



Gl^G 



(70) 



G-^G-J 



and also 



v{q,P) = y{q,P) + 2z{q,P). 
R{P)= [ riq',P)dJ:{q',P) 

Jq' 

X{P)= [ x{q\P)S^q',P) (71) 

Jg' 

Y{P)= I y{q',P)6G{q',P) 

Jq' 

Z{P)= I z{q',P)SG{q',P) 

Jq' 

Note that the condition detM. — from the homogeneous equation [y — 0) gives us back 
fix and 0.2 discussed above. 

In the discussion which follows we look for the bound state referred to above by holding 
Q < Qi < il2 so that it is not necessary here to include the usual scattering ie in the 
denominators. We can invert M. obtaining 



5G^ 



XQrR 



{0^ - 01) 



X 



skv + Bx 



17 



xox 



2 ,2^4ski^+x(]2 yY + 2zZ 



5E 



\sKrR 



+ 



SMX ^^^2 



Bsmx + vQ'^ 



Av + smx 



XX 



Sn 



-A0 

-X(f) 



nx - A{Y + 2Z) 
n{Y + 2Z) - BX 



■{¥ + 2Z) - XVt 



yY + 2zZ 



(72) 



Now we substitute ) in the definitions of tlie quantities R, X, Y and Z ( |7ll) . Because 
r{q, P) is an odd function i.e r{—q, P) = —r{q, P) it is easy to clieck we end up witli a linear 
and homogeneous system for X, Y and Z , that looks like 



W{P,Q).F = 



(73) 



where, omiting the P dependence we have 



w 



W2,\ W2,2 W^2,3 



Y 



(74) 



where the elements of W are given in the appendix B. The system ([73|) will have a non-trivial 
solution if 



detiy(P,(](P)) = 



(75) 



so that we have for each P the corresponding energy ri(P). Numerically the problem reduces 
to calculating determinants of a 3 x 3 matrix. 



A very general property |T8| of the dispersion relation f2(P) can be proven for the par 
ticular case where P = 0. In this case the first line of the matrix ( fflD using ( ^OD and (|64| 



IS 
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A{Q) JgSK{q,0) 



and now using (|39|)-(^0D we have that 



r x\q,0) _X r 1 _ A{0) 
J^SK{q,0) 2 J, JA{q)B{q) A(0)-5(0)- ^ ^ 



So that the first hne of the matrix is zero, making the determinant vanish, showing that we 
will always have a gapless dispersion relation 
independently of the value of A, 

f](0) = 0. (78) 

This zero mode of the RPA equations is the standard Goldstone mode as it's structure 
is associated with the symmetry breaking (indefinite particle number) by the trial wave 
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We note that the dispersion relation that comes from the RPA depends on the total 
momentum defined in ( ^9]) which means that Vt is function of {k — k'), so that it takes into 
account that to obtain an excitation we remove a particle with momentum k' and create 
a particle with momentum k. This result is very different from the Bogoliubov excitation 
where k' is taken to be zero. In other words the Gaussian approximation takes into account 
the effect of the depletion on the excitations. 

B. The Exact Solution 

All the results summarized in this section were obtained by Lieb [|^- [^] who calculated 
the exact ground state energy and also the excitations in terms of independent particle and 
hole excitations. Our purpose is to make a connection with this work which showed that 
interacting bosons in one dimension can be analogous to a fermi gas. We will show how our 
modes correspond to Lieb's particle-hole excitations. 



First of all when A ^ oo it is possible to recover the well-known result |^T| that for an 
infinite coupling constant interacting bosons behave like a system of free fermions, so that 
the ground state energy can be easily computed by 

19 



E 1 1 r^i 

- = -— / e(k)dk (79) 



where Kf is the fermi- momentum. In this case it is trivial to calculate 



1 r^f 

p=— dk (80) 



which gives Kf — np. So that in this particular limit 

E ttV' 



(81) 



AT 2m 3 

We can divide the particle hole excitations of the free fermi gas into two parts. The 
particular hole excitations that correspond to removing a particle from an occupied state to 
just above the fermi level {K^) 

where Kf < k' < Kf and the particle excitations where we remove a particle from the fermi 
level to an unoccupied state 

ejk) = — 83) 

^ 2m 2m ^ ^ 

where k > K{ or k < —K{. To produce a particle hole excitation we must add these 
excitations giving us 

/ 1 1 1\ ^ k Ti k 

eJk,k') = 84 

' ^ 2m 2m ^ ^ 

Using eh and ep one can look at E{k', k) as a function of the total and relative momentum 
P — k — k' and q= {k + k')/2. Because the hole excitations are limited to {—K < k' < K) 
we need to separate two cases, for P > we have 



if P < 2Xf 



if P > 2Ki 



k = Kf 

k' ^Ki-P 

(85) 

k = P-Ki 
k' = -Ki 
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The first case means fixing the momentum of the particle and moving the momentum of the 
hole while in the second one we fix the momentum of the hole at it lowest possible value and 
move the momentum of the particle. So that once we know eh and ep the threshold curve 
defined by the lowest value of E{k', k) for a given P wil be given by 

^ ^ \ e^{Ki - P) + e^{Ki) ioi P < 2K, ^ ^ 

E{P) = \ (86) 

[ eh(-Kf) + ep(P - K,) for P > 2K,. 

Which gives us 

I 27rpP - P2 for P < 2Ki 
E{P) = (87) 
I P2 + 27ipP for P > 2K{ 

Note that the E{P) curve contain two parts. Finite range part P < 2Kf where the con- 
tribution comes basically from the hole excitations. Infinite range part P > 2Kf from the 
particle excitations. The more dilute the system the less finite range part in E{P). 

The generalization for bosons interacting with a finite A was carried out by Lieb and 
Lininger The showed the ground-state can be calculated as 

E 1 



where f{k) is the solution of 



N p J-K -f**''**'* 



with 7 given in (^Tj). The condition 

f^f{k)dk = p (90) 

determines K. For the excitations Lieb defined two different basic interactions which he 
called "particle" and "hole" excitations. To determine these excitation Lieb showed that it 
was sufficient to solve new integral equations. For the particle energy 

^pik) = - + — / pJpip)dp 91 

2m m J-K 
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where Jp{p) could be obtained by solving 



27rJp(p) = 27P 



K 



Mr) 



-K ^^7^ + {p — rY 
and for the hole energy eh(P, A) he had 



dr — 71 + 2 tan ^ 



k — p 

IP 



(92) 



^^ik)=P-^ + -J_^PMp)dp. 



(93) 



where Jh(p) was obtained from 



27rJh(p) = 27P 



K 



Mr) 



-dr + TT — 2 tan ^ 



k' — p 

IP 



(94) 



-K p^^"^ + {p — r)'^ 

In the limit A ^ oo we can see that J(p) and /i (h^ K"^) / (2m) = (^^7r^p)/(2m) 
recovering the free fermion results. Looking at the expressions for the ground state energy 
and the excitations for a given A, Lieb interpreted them as those of a quasi- fermi gas where 
K is an interaction dependent fermi momentum and the distributions factors f{k) and J{k) 
given a special weight for each k. Using this analogy is very reasonable, since it is correct 



in both the A = and A = oo. To obtain the threshold curve E(P) we can use (|86|). The 
difference is that now ep and eh will have different curvatures depending on the interaction. 
Note that in his original work Lieb compared the Bogoliubov scheme with particle and holes 
excitations separately and got very good agreement with the particle excitations. This is 
interesting, as we pointed out earlier, and just tells us that ithe Bogoliubov scheme does 
not contain hole excitations. In general as 7 increases the contribution of the holes for the 
particle-hole excitation energy becomes more and more important and this effect is in part 
described by the Gaussian theory. 



IV. NUMERICAL RESULTS AND CONCLUSIONS 

For the numerical computations we follow Lieb and use the dimensionless coupling con- 
stant 7 and scale all lengths by p and all energies by {h^ p^) / {2m) . In these units we can 
write the ground state energy per particle as 
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^=9{l). (95) 

For the Gaussian static results we solve the non-linear system (^OD-(^OD and determine a 
and h for < 7 < 10. After computing K defined in (|4^) one can get ^'(7) using (|4^) . In 
the Bgoliubov scheme (|46|) leads to 



^b(7) = 7 



(96) 



Finally the exact result for (7(7) was obtained by Lieb solving (|88|)-(]90|). The results as a 
function of 7 can be seen in Fig. 1. Where one can basically see the very good agreement of 
both Gaussian and Bogoliubov) for low 7 (7 < 1) with the difference that by construction 
the Gaussian result It is always an upper bound and the Bogoliubov energy is below the 
exact result. For higher values of 7 the Bogoliubov results, as expected, collapses (7 ~ 5) 
while the Gaussian theory still gives a result. In the hard core limit (7 — 00) the Gaussian 
approximation fails to go to the finite free fermi energy of the exact solution. 

To obtain the "exact" excitation energies we used Lieb's results for eh and obtained 
by solving (pTD-(p^). Using (|86| ) we obtain the excitation curves E{P) as a function of P as 



shown in Fig. 2 for 7 = 0, 0.787, 4.527 and 00. Note that for P < 2K the contribution from 
the hole excitations produce quite interesting dispersion curves. 

To illustrate the RPA results we plot in Fig. 3, Qi, Q2 and f2 as a function of P for 
7 = 7.551. The relative momentum variable q will gives us a continuum that corresponds 
to a scattering region and one can see clearly how Qi gets pushed down by the interactionp 
producing the gapless fl. We compare the bound state dispersion relation Q with the 
Bogoliubov one and with the exact results for 7 = 0.787 (Fig. 4) and 7 = 4.527 (Fig. 5). 
We note that the improvement of the RPA compared to the Bogoliubov result increases with 
7. When 7 goes to infinity both the RPA and the Bogoliubov schemes fail. 

The results of this paper can be summarized in three points 

1) We have seen that the Gaussian variational method takes into account the self- 
interactions of the particles out of the condensate. This fact starts improving the ground 
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state energy results when compared with Bogohubov theory once 7 increases (7 > 2). It is 
very difficult however for the Gaussian variational results to get close to the exact for high 7 
values. In the one dimensional case when 7 — > 00 the Gaussian ground state energy diverges 
while the exact goes to the finite fermi energy. 

2) For the particle hole excitations the Gaussian results take into account that when we 
remove a particle from an occupied state in order to make an excitation, this particle can be 
out of the condensate which does not happen in the Bogoliubov scheme. Again this seems 
an improved description for intermediate values of 7. As 7 increases the depletion starts to 
grow and in this particular case goes to a fermi sea when A ^ 00 which can not be described 
by our methods since they do not deal with very short range correlations. 

3) Another important point that is somehow related to item 2 is that to solve the RPA 
equation we had decoupled harmonic oscillators which , in the quasi-boson picture means 
that the RPA takes into account the interaction between the quasi-bosons generating new 
modes that are mixtures of one (50, Sir) and two quasi-bosons {5G, 5E). In the Bogoliubov 
scheme we always have free quasi-bosons when their interaction starts to be relevant as 7 
increases. 

We can conclude that the Gaussian variational method can describe systems where the 
depletions can not be neglected (when dilute theories break down) but because of the absence 
of short range correlation will require corrections for highly depleted systems (like Helium 
4). 
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APPENDIX A: 



If we carefully examine the expression (|T9|) we see that the functional integral that must 
be computed involve moments of a Gaussian. Namely, 

J (P(/.')^l0',t]0'(xi) . . .0'(x4)v&[0',t] (Al) 

where "^[(f)'] is our normalized trial wave functional given by (|^). This functional integrals 
can be computed easily if we include a source term [E^ in the normalization integral i.e. 



Using the expression for \E'[0',t] we complete the squares in the exponential getting 

The source term allows us to compute the functional integral of any moment of the Gaussian 



J=0 



For instance we can calculate 



t]</.'(x)</.'(y)y (x - y)0'(x)0'(y)M/[0', t] = 40(x, t)<P{y, t)V{^ - y)G((x, y, t) 
+</)2(x, t)y (x - y)G'((y, y, t) + ^^(y, t)V{^ - y)G'((x, x, t) 

2G'(x, y, t) V(x - y)G'(x, y, t) + ^(x, x, t)\/(x - y, t)G{y, y, t) (A5) 



Note that the last two terms of ( |A5|) show, as expected, the mean field factorization. 



APPENDIX B: 



After some algebra one gets the elements of the matrix W which appears in equation 
[[f5|). For calculating the dispersion relation, Q will be below the lowest pole appearing in 
(B1)-(B9). 
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W^3,2 = -A^0^ (02_J^2)(^2_J^2) - ^ (02 - ^i) ^^^^ 
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Figure Captions 



FIG.l The ground state energy per particle can be written as E/N = p'^g{'y). The curves 
give (7(7) as a function of 7 for the Gaussian , the Bogohubov and the exact resuh. 

FIG.2 The exact threshold curves for the particle-hole excitation as a function of the 
total particle - hole momentum P = k — k' as obtained using (RBf) 



FIG. 3 The free quasi-boson energy VLi[P) given by (|63|). The two free quasi-boson 
energy Vt2{P.,q) given by (|64D and the new gapless ^{P) mode energy obtained with the 
RPA calculation by solving (^). Note that, as expected, the gapless property of ^{P) 
comes from the Goldstone mode associated with the symmetry breaking 

FIG. 4 Comparison between the three different particle-hole excitations curves: Bogoli- 
ubov ^\{P) given by (|66D, the RPA VL{P) by solving ( ffSD and the exact using (pB)). These 
are given as a function of the total particle-hole momentum P for 7 = 0.787 . 

FIG. 5 Comparison between the three particle-hole excitations curves: i^i(-P) given by 
(|66|), ^{P) by solving ( fTSD and the exact using ( p6D as a function of the total particle-hole 
momentum P for 7 = 5.527. 
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